Emergence of foams from the breakdown of the phase field crystal model 
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The phase field crystal (PFC) model captures the elastic and topological properties of crystals with 
a single scalar field at small undercooling. At large undercooling, new foam-like behavior emerges. 
We characterize this foam phase of the PFC equation and propose a modified PFC equation that 
may be used for the simulation of foam dynamics. This minimal model reproduces von Neumann's 
rule for two-dimensional dry foams, and Lifshitz-Slyozov coarsening for wet foams. We also measure 
the coordination number distribution and find that its second moment is larger than previously- 
reported experimental and theoretical studies of soap froths, a finding that we attribute to the 
wetness of the foam increasing with time. 
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Computational modeling of large-scale systems usu- 
ally involves either detailed molecular dynamics simula- 
tions, in which every particle must be tracked, or a highly 
coarse-grained model in which the underlying symmetries 
are known and are used to derive a continuum model 
from the underlying physics. Molecular dynamics is lim- 
ited to small systems and/or to very short times, and 
coarse-grained models tend to fail at places where the 
underlying symmetries are broken, such as at defects and 
dislocations. The phase field crystal (PFC) model pQ is 
an intermediate approach with the advantage of diffusive 
time scale, but with atomic scale resolution of molecu- 
lar dynamics. The PFC model can be used to capture 
the dynamics of defects and dislocations in large crystals 
[2j |3] , the dynamics of grain interactions [4 , molecular 
dynamics of vacancies [5] and even nonlinear elasticity [6]. 
In addition, it is amenable to methods such as coarse- 
graining and adaptive mesh refinement [7]. 

The dimensionless phase field crystal equation [1 

^ = V 2 [(V 2 + l)^ + ^ + ^ 3 ] (1) 

is a density functional theory [2], best thought of as 
arising from phenomenological and symmetry consider- 
ations. In particular, this is the simplest class of models 
appropriate for systems whose dynamics is governed by 
minimizing departures from periodicity [8], as opposed to 
the situation in other materials processes, such as spin- 
odal decomposition, where the dynamics minimizes de- 
partures from spatial uniformity. In principle this ap- 
proach only holds for small values of the undercooling 
a = — r, beyond which the strong nonlinearity may in 
principle overwhelm the crystal's symmetries, leading to 
such artifacts as merger or dissolution of the 'atoms' of 
the model. 

In this Rapid Communication, we explore an interest- 
ing feature of the phase field crystal model in the limit 
of large undercooling, which we show turns out to facili- 
tate a simple, scalar and minimal model of foams. In the 



appropriate density regimes, the behavior of the phase 
field crystal equation at large undercoolings is that the 
atoms of the crystal lattice begin to merge. However, the 
interstitial spaces between the atoms are preserved and 
become line solitons. The result is a coarsening foam-like 
structure. We analyze the process by which this insta- 
bility in the equation of motion occurs, and use this un- 
derstanding to propose a minimal, modified PFC model 
that is capable of describing quantitatively both wet and 
dry foams at the level of a continuum, scalar theory that 
is computationally efficient. 

Equilibrium Phase Diagram:- The phase diagram 
of the PFC model has been computed for small 
undercoolings [1. We shall use the same methods to con- 
struct the phase diagram at larger values of the under- 
cooling in order to see what may be found there. First, 
however, we will convert the PFC equation and PFC en- 
ergy to a nondimensionalized form with respect to the 
equilibrium liquid (constant phase) density. The energy 
minimizing density for the constant phase is ijj = ±y/a, 
so we will introduce a nondimensional order parameter 
(j> = ip/y/a. This gives us the following free energy: 

F = f ^(V 2 + 1)V + a (-^ + J/) dV (2) 

V 

and corresponding equation of motion: 

^ = V 2 ((V 2 + 1) 2 + «(-</> + 3 )) (3) 

We then construct the phase diagram at large a by 
calculating the energy minima of the one-mode approxi- 
mations for the constant phases L± : (j) = </>o, where the 
subscript ± refers to the sign of the average density 0o; 
the striped phase S : (j) = c^o + As cos(x), where As = 
\A*(4 — 3(/>q)/3; and the two triangular lattices A± : <p = 
0o + (^a±^a)(cos(/ci -r) + cos(A:2 -r) + cos(/c3 -r)), where 
the kj are the lattice vectors of the regular triangular lat- 
tice, A/\ = — acf)o/5 and = ^/a(15 — 360q)/15. 
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Even though the one-mode approximation is not accu- 
rate at these large values of a, we use this approximation 
to understand heuristically the behavior observed, thus 
motivating our form for the modified phase field crystal 
model introduced below. Our main results, for the mod- 
ified phase field crystal model, are fully time-dependent 
and independent of this approximation. 
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FIG. 1: (Color online). Phase diagram associated with the 
free energy in Eq. Invariant reactions occur at a = 2.56 
and a = 5.45. No further topological changes occur for a > 8. 

The equilibrium phase diagram is obtained by substi- 
tuting the various ansatz into Eqn. (|3| and performing 
a common tangent construction. The result is shown in 
Fig [I] We find a series of previously undiscovered invari- 
ant points as a increases. At a = 2.56, the coexistence 
region between the striped and triangular phases disap- 
pears, giving rise to a coexistence between stripes and 
one of the liquid phases. Then, at a = 5.45 the striped 
phase vanishes, leaving only a region of immiscible liquid- 
liquid coexistence, and the regions comprised entirely of 
the liquid phases. 

This can be understood by considering the influence 
of the two terms in the free energy. At small a, the 
wavelength selection term that was added to construct 
the PFC equation is dominant. As a increases, the local 
term that drives <\> to have the values of ±1 becomes 
dominant, so that eventually the wavelength selection 
term (1 + V 2 ) 2 )(/> is a small perturbation and the equi- 
librium phases are determined by the average value of 
0o/2 — </>q/4, which is minimized by the constant phases. 

Dynamics:- The phase diagram we have just computed 
represents the behavior of this system in its final equilib- 
rium state. However, the approach toward that equilib- 
rium changes as a increases. Starting from a hexagonal 
lattice, the system attempts to coarsen into two regions, 
one composed of L + , the other composed of L_. How- 
ever, because of the small residual wavelength selection, 
there is a finite energy cost to removing spatially pat- 
terned structures. The consequence of this is that the 



PFC 'atoms' are dynamically conserved by this infinites- 
imal energy barrier, which halts the coarsening dynamics. 
If the average density of the system is shifted, e.g., by de- 
liberately reducing </>o over time, this can provide enough 
energy to destabilize the PFC atoms, but stripe-like cell 
walls will still remain stable for a longer period as they 
can move perpendicular to their length to create large 
bubbles of one of the liquids. 




FIG. 2: a. Final state obtained from Eqn. (|3| with a — 16 
and 0o = —0.41 after coarsening. Residual atoms coexist with 
large foam cell regions, b. Simulation using the modified PFC 
equation for a — 20, quenched with 0o = 0.4 initially, and 
then drained to 0o = —0.6. The result is a dry foam structure 
with no residual atoms, c. Simulation of the modified PFC 
equation for a = 20 and 0o = 0.25, allowed to coarsen without 
draining. The result is a wet foam structure with circular 
bubbles. 

The end result of this is that a foam-like state appears 
(FigJ^i), which coarsens to a stationary point at which 
there are no more free adjustments that can be made 
to approach the equilibrium state. Thus, even a small 
perturbative addition of wavelength selection is sufficient 
to stabilize a foam-like structure. 

Although this behavior is interesting from the point of 
view of understanding the properties of the phase field 
crystal model, it does not behave like a physical foam. 
The arrested coarsening and inability to get rid of resid- 
ual atoms prevent this from being used as a model for 
studying coarsening foams. We now alter the free energy 
so as to destabilize the residual atoms while retaining the 
overall cellular structure, so that we may recover physical 
foam dynamics. 

In effect what we wish to do is to weaken the wave- 
length selection while encouraging k = structures. A 
stripe has one k — direction and one direction with the 
selected wavelength, whereas an atom has two directions 
with the selected wavelength. If the selected wavelength 
sits at a local energy minimum, but the global minimum 
is at k = 0, then this should help remove residual atoms 
in favor of bubble interfaces. 

The PFC wavelength selection operator (V 2 + l) 2 cor- 
responds in /c-space to 

E(k) = (l-k 2 ) 2 (4) 

We modify this by introducing two parameters ko and b 

E mod (k) = (k/k ) 2 (l - (k/ko) 2 ) 2 + b(k/k ) 2 (5) 
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and then choose ko so as to retain minima at k = ±1: 
k = [3/(2 + vT^S)] 1 / 2 . The modified form of the free 
energy in real space becomes 

wKH4 v2+i ) 2 -H* 

v v 7 

+ a (-\^ 2 + \^) dV ( 6 ) 

We can then vary b to control the relative depths of the 
minima at k = ±1. The most extreme effect is achieved 
at b = —1/3, at which the minima at k = ±1 disappear. 
We use this value of b throughout the rest of this paper. 

The dynamics of the modified PFC equation give rise 
to the structures shown in Figs. [2]o and If the sys- 
tem is quenched from a disordered state, one recovers 
Lifshitz-Slyozov coarsening dynamics [9] corresponding 
to the physical behavior of a very wet foam. On the 
other hand, if the system is drained from a state with pos- 
itive average density to one with negative average den- 
sity, polygonal cell walls form and coarsen by topological 
rearrangements. 

Results:- We measured the coarsening dynamics and 
statistics of the eventual scaling state of the modified 
PFC equation in order to compare with physical foams. 
In order to measure the coarsening we are interested in 
the quantity (r)(t), the average bubble radius. We may 
easily count the total number of bubbles and so determine 
the average bubble area as (A)(t) = A tota i/N(t). If we 
assume that deviations from circular geometry are small, 
then we may approximate (r)(t) ~ \J (A)(t). 
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FIG. 3: Coarsening behavior of wet and dry foams in the 
modified PFC equation with a = 20, compared with the cor- 
responding theoretical coarsening laws. 

For a wet foam, we expect the bubble growth to pro- 
ceed as diffusive grain growth, so that (r)(t) cx t 1 / 3 . As 
an example, we simulate a 1024x1024 system with a = 20 
and quench into an average density </>o = 0.25. We ob- 



serve a scaling (r) ex t 0,34 (Fig[3j filled circles), consistent 
with the predicted growth law. 

For a dry foam in two dimensions, von Neumann's law 
for bubble area growth [10] is 



where k is an effective diffusivity, and n is the coordina- 
tion number of the bubble. This implies that the average 
bubble radius should scale as (r)(t) cx t 1 / 2 . To measure 
the coarsening of a dry foam, we quench from an aver- 
age density </>o = 0.3 and drain slowly to a density of 
0o = —0.4. We start measuring the coarsening dynamics 
after we have stopped draining. However, our reference 
point for t = is still the point of the quench. We observe 
a scaling r cx t 0,43 , significantly slower than the predic- 
tions of von Neumann's law (Fig|3j open circles). At late 
times in the simulation, the bubble interfaces start to be- 
come wet. If we restrict our analysis to times shortly after 
coarsening begins, the resulting fit exponent increases to 
r cx t 0A7 . This indicates that we are probably seeing the 
effect of not having a fully dry foam at any point in time. 

As two-dimensional froths coarsen, they are expected 
to reach a self-similar scaling state in which the normal- 
ized moments of the distribution of areas and of coordi- 
nation number are expected to become time independent. 
The second moment of the coordination number distri- 
bution has been used as a probe of this scaling state. 
Glazier and Weaire [IT] predict that the coordination 
number distribution should eventually reach a universal 
limiting scaling state with a second moment ji2 = 1-4, 
with a strongly non-monotonic transient behavior. 

We measure the coordination number distribution of 
the dry foam by a watershed algorithm. [12] We super- 
impose a grid over the system and fill each bubble with a 
unique integer. We then allow these integers to propagate 
to neighboring unoccupied cells, and then find the total 
number of different integers that a given bubble comes 
into contact with. We observe a similar time-dependence 
of the second moment of our distribution, but the mea- 
sured limiting value H2 = 2.31±0.01 is significantly differ- 
ent from 1.4 (see Fig|4|. In other models and simulations 
of coarsening 2D foams, the observed limiting value of \±2 
has been variously reported as 1.9 [13;, and 1.7 [14 , and 
1.5 [15], and experiments have measured from 1.4 [16] in 
soap froths to values as small as 0.14 and 0.22 in magnetic 
froths HZJ. 

We conclude that the second moment is either not a 
universal quantity, or else that there are strong transients 
which make any universal scaling regime difficult to ob- 
serve in practice. We note that our foam is not com- 
pletely dry, and that the absence of a drainage mech- 
anism implies that the foams becomes wetter the more 
they coarsen. Whether real foams undergo a correspond- 
ing change of regime is not clear. 

We have found that in the limit of large undercooling, 
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FIG. 4: Evolution of the second moment of the bubble co- 
ordination number in a dry foam realized by the modified 
PFC equation. The bubble distribution broadens during the 
transition to power-law coarsening, after which it reduces to 
a steady value of 2.31 =b 0.01. The inset shows a particular 
distribution at t = 1800. 



the phase field crystal equation produces topologically 
stabilized foams. These foams are a consequence of the 
residual wavelength selection, which acts as a singular 
perturbation that prevents the destruction of cell walls 
and forces the coarsening to take place via topological re- 
arrangements. The fact that this behavior emerges from 
a minimal model such as the phase field crystal equation 
suggests a general mechanism by which foams may occur 
in natural systems. 

The ingredients of the phase field crystal equation are a 
driving force towards certain equilibrium densities, some 
sort of spatial relaxation (diffusion, viscosity, and the 
like), and some competing source of wavelength selection. 
Such a system, in the limit where the wavelength selec- 
tion is weak compared to the other forces, would likely 
give rise to a foam. This mathematical structure can be 
seen in models of magnetic froths [17] . Type-I supercon- 
ductors [18] , and in models of polygonal cells found in 
melting snow [T9] , 

However, the foams produced in this limit of the PFC 
equation have unphysical properties. The wavelength se- 
lection also preserves atoms: bubbles whose diameter is 
comparable to the width of the cell walls. This can be 
ameliorated by modifying the PFC free energy to pe- 
nalizes atoms while encouraging stripes. As a result, the 
coarsening dynamics of real foams are recovered although 
the distribution of bubbles in the resultant scaling state 



appears to be somewhat different. 

If these differences are understood, the modified PFC 
equation may be a useful tool in modelling foams, as it 
is simpler than other existing methods such as Q-Potts 
models [20] and minimal surface evolution [21], and is 
fairly easy to simulate, having only one field and a spatial 
structure that is easily treated with spectral methods. 
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